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We present the new nCTEQlB set of nuclear parton distribution functions (nPDFs) with uncer¬ 
tainties. This ht extends the CTEQ proton PDFs to include the nuclear dependence using data on 
nuclei all the way up to ^°®Pb. The uncertainties are determined using the Hessian method with 
an optimal rescaling of the eigenvectors to accurately represent the uncertainties for the chosen 
tolerance criteria. In addition to the Deep Inelastic Scattering (DIS) and Drell-Yan (DY) processes, 
we also include inclusive pion production data from RHIC to help constrain the nuclear gluon PDF. 
Furthermore, we investigate the correlation of the data sets with specihc nPDF flavor components, 
and asses the impact of individual experiments. We also provide comparisons of the nCTEQlB set 
with recent fits from other groups. 
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I. INTRODUCTION 

In the last thirty years, an impressive array of dis¬ 
coveries in particle physics has come from high energy 
hadron experiments. These discoveries, along with many 
other key measurements, rely on our understanding of 
nucleon structure. A nucleon can be described using the 
language of parton distribution functions (PDFs) which 
is based on QCD factorization theorems [TH3]. PDFs 
are determined in global analyses of a variety of different 
hard scattering processes such as deep inelastic scatter¬ 
ing (DIS), Drell-Yan (DY) lepton pair production, vec¬ 
tor boson production and the inclusive jet production. 
The backbone of any global analysis are the very precise 
DIS structure function data from HERA which cover a 
wide kinematic range in {x, Q^)- Several global analyses, 
based on an ever growing set of precise experimental data 
and on next-to-next-to-leading order (NNLO) theoretical 
predictions, are regularly updated and maintained I2HS]- 
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Over the years, a series of global analysis studies have 
been performed within a single framework, or comparing 
different frameworks. For example, detailed studies of 
PDF uncertainties have been compared using Hessian, 
Lagrangian and Monte Carlo methods. Furthermore, the 
precision of experimental data and theoretical predictions 
in the proton case allows one to perform studies of smaller 
effects such as the difference between the treatment of 
heavy quarks in different analyses or the exact treatment 
of target-mass corrections and higher twist effects. As a 
consequence, the nucleon structure is quite well known 
over a wide kinematic range. 

Similarly, the theoretical description of hard scattering 
processes in lepton-nucleus and proton-nucleus reactions 
requires the knowledge of parton distribution functions 
inside nuclei characterized by the atomic number A and 
the charge Z. It has been known since the discovery 
of the EMC effect [10] more than 30 years ago that the 
nucleus cannot be considered as an ensemble of Z free 
protons and (A — Z) free neutrons. Consequently, the 
nuclear PDFs (nPDFs) will differ from the naive addi¬ 
tive combination of free proton and neutron PDFs. As 
in the proton case, nuclear PDFs have been determined in 
the literature by global fits to experimental data for hard 
scale processes including deep inelastic scattering on nu¬ 
clei and nuclear collision experiments [IIHIll- However, 
compared to the proton our knowledge of nuclear PDFs 
is much less advanced. There are several reasons. 

On the theoretical side, the description of nuclear in¬ 
duced hard processes is more challenging due to the com¬ 
plex nuclear environment. Still, all global nuclear PDF 
analyses rely on the assertion that the QCD factoriza¬ 
tion theorems remain valid for (.A and pA hard scattering 
processes, see e.g. miiig. In fact, it is only in this con¬ 
text that the universal parton distributions {f^{x,Q)) 
are defined; they are given as matrix elements of the 
same local twist-2 operators as in the proton case but on 
nuclear states. The nuclear PDFs then account for nu¬ 
clear effects (in particular EMC suppression, shadowing, 
anti-shadowing) at the twist-2 level in a universal manner 
and the entire formalism becomes predictive. However, 
higher twist contributions are expected to be enhanced in 
a nucleus (oc A^/^) [121111]. Here, final state re-scattering 
corrections due to the propagation of the outgoing par- 
tons through the nuclear medium, which are higher twist, 
should be power suppressed but may be substantial and 
so must be either included in the analysis or eliminated 
by suitable kinematic cutsQ In addition, other effects like 
a different propagation of the hadronic fluctuations of the 
exchange bosons in the nuclear mediumj^ gluon satura- 


^ Needless to say that the final state interactions do not concern 
the fully inclusive DIS structure functions but may be relevant 
for less inclusive observables (single pion production, di-muon 
production in vA DIS, ...). On the other hand, power sup¬ 
pressed initial state interactions are expected to be numerically 
small. 

^ There could be modifications of charged current neutrino scat- 


tion, and deviations from DGLAP evolution at small x 
may play a more prominent role in the nuclear case, see 
e.g. [IHldS] and references therein. 

Ultimately, the validity of the twist-2 factorization for¬ 
malism will be tested phenomenologically by how well 
our approach based on the factorization assumption de¬ 
scribes the data. The existing global analyses generally 
lead to a good description of the data confirming this pic¬ 
ture; however, it may be challenged by future precision 
data from the LHC and an Electron-Ion Collider (EIC) 
covering an extended kinematic plane. It is notable that 
tensions between vA DIS data and lA DIS data have 
been reported [23 HI] which might be due to higher twist 
contributions, or indicate a breaking of twist-2 factoriza¬ 
tion. These tensions largely disappear if the correlations 
between the NuTeV data points are discarded [22] • 

The other reason why nuclear PDFs lag behind the 
proton analyses can be traced back to the lack of pre¬ 
cise experimental data. For example the constraints on 
the nPDFs for any single nucleus are (so far) too scarce, 
so that experimental data from scattering on multiple 
nuclei must be considered. Since the nuclear effects are 
clearly dependent on the number of nucleons, this re¬ 
quires modeling of the non-trivial nuclear A dependence 
of the parton distributions. Even after combining the 
data sets for different nuclei, the precision of the nuclear 
PDFs is not yet comparable to the proton PDFs where 
quark distributions for most flavors together with the 
gluon distribution are reliably determined over a broad 
kinematic range, due to the smaller number and hence 
smaller kinematic coverage of the current relevant nu¬ 
clear data. As a consequence, the nuclear PDFs in every 
analysis have large uncertainties as the parton distribu¬ 
tions are not fully constrained by the available data. The 
nuclear PDFs still largely depend on assumptions inher¬ 
ent in every analysis. The dependence on assumptions, 
such as for example the parameterization form, leads to 
predictions where different analyses differ by more than 
the estimated uncertainties. It follows therefore that in 
order to assess the true uncertainty, all available results 
and their uncertainties should be considered and com¬ 
bined. 

In this paper we present a new analysis of nuclear 
PDFs in the CTEQ global PDF fitting framework. We 
use theoretical predictions at the next-to-leading order to 
fit all available data from charged lepton DIS and Drell- 
Yan di-lepton production as in our previous analysis [la¬ 
in addition, we have added inclusive pion production 
data from RHIC and have performed a careful analysis of 
the uncertainties using the Hessian method. Our frame¬ 
work differs considerably from other global analyses of 
nuclear PDFs which we compare our results with. 

The remainder of this paper is organized as follows. In 


tering that are different than those for neutral current charged 
lepton scattering for instance due to the exchange of a charged 
massive vector boson El- 
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Section [T^ we introduce in detail the framework includ¬ 
ing the parametrization of the nPDFs at the input scale 
together with a review of the Hessian method which we 
use to estimate the uncertainties on the nPDFs. In Sec¬ 
tion III we review the experimental data included in the 
fit. In Section |IV) we present the results of our fit, com¬ 
pare with recent results from the literature, and examine 
the correlations between individual PDF flavors and the 
various experiments. Finally, in SectionjVj we summarize 
the obtained results. Additionally we include two appen¬ 
dices. In Appendix [A| we provide details on the Hessian 
rescaling method, and in Appendix we comment on 
the usage and availability of our nPDFs. 


II. THE NPDF FRAMEWORK 

In this section we describe in detail the framework 
of the nCTEQ global analysis. For the purpose of fit¬ 
ting nuclear parton distributions we will parameterize 
the PDFs of a proton bound in a nucleus 
A, then construct the full distributions of partons in the 
nucleus using isospin symmetry, and in the end perform a 
fit just like in the case of the free proton. Indeed, isospin 
symmetry is used to construct the PDFs of a bound neu¬ 
tron, Q)i from those of the proton by exchanging 

up- and down-quark distributions. Afterwards the par- 
ton distributions of the nucleus are constructed as: 

#’^)(x,Q) = + ^^/”/^(x,Q), (2.1) 

where Z is number of protons and A number of protons 
and neutrons in the nucleus 0 

The theoretical calculations in our global analysis 
make use of parton distributions of a particular nucleus 
determine the DIS structure functions, Drell- 
Yan cross sections or the cross section for an inclusive 
pion production: 

Q') = E Q') ® g'), (2.2) 

d<XAB^llX = E ® ® 

(2.3) 

= E ® ® ® Dl, (2.4) 

ijk 


^ Note that the PDFs of the nucleus, are the ob¬ 

jects of interest which are constrained by the experimental data, 
whereas the Q) and Q) are just effective quanti¬ 

ties used internally to decompose the nuclear PDFs. They should 
not be interpreted literally as matrix elements of local operators 
where the free nucleon states have been replaced by bound nu¬ 
cleon states in a nuclear medium since they also include effects 
from multi-nucleon states. The notion of “effective bound nu¬ 
cleon PDFs” is also used in the literature discussing the factor¬ 
ization in the case of pA interactions US]. 


where 0 stands for a convolution integral over the mo¬ 
mentum fraction. The DIS structure functions calcu¬ 
lations are carried out using the ACOT variable flavor 
number scheme 0 [13H15] at next-to-leading order in 
QCD [5S] We take into account only the dominant tar¬ 
get mass effects which are included in the structure func¬ 
tion expressions in the ACOT scheme [53]. Full treat¬ 
ment of the target mass corrections |30] is not necessary 
in our analysis because they are relevant mostly at large 
X and low a region of phase-space which we exclude 
by kinematic cuts. Moreover, the target mass corrections 
are expected to be of lesser importance in the ratios of 
structure functions. 

In all theory calculations we identify the renormaliza¬ 
tion and factorization scales: p. = = fj.p. The scale 

is set differently for different processes: in deep-inelastic 
scattering it is set to the virtuality of the exchanged vec¬ 
tor boson in Drell-Yan production processes it 

is set to the invariant mass of the produced lepton pair 
and in inclusive pion production the common 
scale is set equal to the final state fragmentation scale as 
p. = fi'p = 0.5pT where px is the transverse momentum of 
the produced 7r°. To speed-up the evaluation of next-to- 
leading order cross sections in the fit, we have the ability 
to use K-factors; however for the final fitting the full 
NLO calculations are used. In the case of inclusive pion 
production, we use the results of Ref. isiiisa and speed 
up the calculation by using pre-computed grids already 
including convolutions with one PDF and fragmentation 
function and leaving only one convolution (with the nu¬ 
clear PDFs) to be calculated during the fitting procedure. 


A. Parameterization 

The starting point of any determination of parton dis¬ 
tribution functions is the parameterization of individual 
distributions at the input scale Qo- The parameteriza¬ 
tion of the presented nCTEQ nuclear PDFs is the same as 
in our previous analyses [iliiilisg. It mimics the pa¬ 
rameterization used in the free proton CTEQ fits |33H53] , 
and takes the following form: 

x/f/^(x, go) = CO x^Hl - 

for i = Uy,dy,g,u + d, s + s,s — s, 

= Co x'^Hl - x)=" + (l-bC3x)(l-x)‘=F 

u{x,Qo) 

(2.5) 

The input scale is chosen to be the same as for the free 
proton fits |34l|36|, namely Qo = 1.3 GeV. 


^ For recent extensions of the ACOT scheme to higher orders, re- 
quired for global analyses at next-to-next-to-leading order, see 
[27[I28| : the massless limits have been validated with the help of 
QCDNUM.[29] 
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However, this parameterization needs to be appro¬ 
priately modified to accommodate the additional nu¬ 
clear degrees of freedom. As in other available nuclear 
PDFs [TTJ [TH HU , nuclear targets are characterized only 
by their atomic mass number A. However, in contrast to 
those nPDFs where the nuclear effects are added on top 
of the free proton PDFs in form of ratios, in our analysis 
we introduce the additional A dependence directly to the 


c-coefhcients of Eq. (2.5): 


Cfc —t Ck(A) = Ck,o -\- Cfcp (l — A , 

k = {1,...,5}. 


( 2 . 6 ) 


This parameterization is designed in such a way that for 
A = 1 one recovers the underlying PDFs of a free proton. 
The free proton PDFs are described by the coefficients 
Cfc^o which in our analysis are fixed to values of the fit 
of Ref. [31] which is close to CTEQ6.1 |3S| but has the 
advantage of having minimal influence from nuclear data. 

Although in principle this framework can be used to 
determine the strange quark content of the bound nu¬ 
cleon, there is not sufficient data available to reliably do 
that. Therefore we assume that at the initial scale Qo 

sP/^{x,Qo) = sP/^{x,Qo) = + , (2.7) 

where k{A) is an A-dependent normalization factor pa¬ 
rameterized as k{A) = -I- Cg()'*(l — A“'^o^ ® 

The normalization coefficients cq in Eq. ( |2.5| ) are differ¬ 
ent than the other parameters. They are also dependent 
on the atomic number but not all of them are free param¬ 
eters that can be fitted. Most of them are constrained by 
sum rules. The normalization coefficients for the valence 
quark PDFs are constrained for each atomic number A 
by requiring that they obey the number sum rules 

[ dx Qo) = 2 , [ dx fjf'ix, Qo) = 1 ■ 

( 2 . 8 ) 

The remaining normalization coefficients are constrained 
by the momentum sum rule 


of s -I- s to be determined durin g th e glob al fit together 
with the parameters from Eqs. (|2.5|) and (2.6). The A- 


dependent momentum fraction of gluon is parametrized 
as 


dxxgPl^ix, Qo) = Mg exp 


c5,o + ci 


(l- a-<2^ 


(2A0) 

which modifies the momentum fraction of the gluon in a 
free proton (described by coefficients Mg and Cg g). 

The momentum fraction of the s -|- s combination is 
then given by 


J dx x(^sP^"^{x,Qo) 

K 

(2 -|- k) 


aP/A 


(a^,Qo)) 


( 2 . 11 ) 


i'-L ""S 




nS + S 
-0,0 


s-l-s ( 

0,1 


I - A 


where the sum runs through i = Uv,dy, g. The remaining 
normalization parameters are taken care of by the mo¬ 
mentum sum rule and do not introduce additional free 
parameters. 

The parameterization of Eq. ( |2.5| ) together with the 
whole nCTEQ nuclear PDF framework has been designed 
in analogy to the free proton PDFs where parton momen¬ 
tum X is restricted to be in the range (0,1). However, in 
the nuclear case, x represents the parton fractional mo¬ 
mentum with respect to the average momentum carried 
by a nucleon. Since a particular nucleon can have a mo¬ 
mentum bigger than an average nucleon, x can extend 
up to A in a nucleus with an atomic number A. If one 
were to take this into account, one would have to modify 
the sum rules in Eqs. (2.8) and (2.9) together with the 
DGLAP evolution. However, the structure functions at 
a: > 1 fall off rapidly and the contribution to the mo¬ 
ments of the structure functions from the region of a: > 1 
is very small ISZlISHj. Therefore, all currently available 
nuclear PDFs have been obtained neglecting the a: > 1 
region and we follow the same pathj^ 


B. Finding the optimal PDFs 


[ dxY^ xfj^{x, Qo) = 1, 


(2.9) 


which however can only determine one of them. The rest 
of the normalization parameters are either considered as 
free parameters in the fit or are fixed using additional 
assumptions to simplify the analysis (e.g. like Eq. (2.7|). 
We choose to introduce free parameters for the momen¬ 
tum fraction of the gluon and for the momentum fraction 


The fitting procedure used to find PDFs that describe 
the considered data best is based on minimizing the ap¬ 
propriate function, as described in [33]. The simplest 
definition of the x^ function for n experiments is 




[A — Ti{{aj})]'^ 


( 2 . 12 ) 


where Di are the measured experimental values, Ti are 
the corresponding theoretical predictions and erf are the 


^ This is a straightforward generalization of the approach employed 
in the underlying proton analysis which also assumes that at 
the initial scale Qo the strange quark PDFs are constrained by 
s = s = d). 


® In fact the first next-to-leading order nuclear PDF analysis m 
used a framework which at least in principle allows to accommo¬ 
date the case of a: > 1. 
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systematic and statistical experimental errors added in 
quadrature. The parameters {aj} are a set of free pa¬ 
rameters which define the PDFs at the input scale (see 
Eq. (2.5)) and are varied in order to find the minimum 
of the function. 

This simple definition, with slight modifications al¬ 
lowing for the inclusion of overall changes to data normal¬ 
ization, is used by most of the groups performing nuclear 
global analyses. However, in the current analysis, as in 
the previous nCTEQ fits [laEniEiiiis] this simple dehni- 
tion is modified to account for correlations in the exper¬ 
imental uncertainties. We follow here the prescription 
suggested in Ref. [35]. The total for n experiments 
with parameters {a^} is defined to be 


X^{{aj}) = xl{{aj}) , (2.13) 


where is the weight for experiment n; for our fits all 
weights are set to 1. The Xn ^ contribution from one 
individual experiment n, and this is given by 




^ [A - Ti({aj})]'^ 


A Bk' 

k,k' 


(2.14) 

where i runs over data points and k, k' run over sources of 
the correlated uncertainties. For each experimental data 
point we sum the statistical error ai together with the 
uncorrelated systematic error Ui in quadrature to obtain 
= erf -I- uf. The components of the correlated uncer¬ 
tainties are given by [3S| 


Bk{{aj}) = 


P^k [A-T,({a,})] 


z/cfc' — 5kk' + 'y ^ 


PikP: 


ik' 


(2.15) 


where j3ik are the sources of correlated systematic errors. 

We stress that in this procedure only the experimental 
uncertainties are accounted for; all theoretical and model 
uncertainties (e.g. missing higher order corrections, pa¬ 
rameterization choice, etc.) are not taken into account. 

Having defined the appropriate x^ function it needs 
to be minimized with respect to the fitting parameters 
{aj} that define the bound proton PDFs at the ini¬ 
tial scale Qq. We perform the minimization using the 
pyMinuit package [ID] which is a python interface to 
“SEAL-Minuit” [IT] — a C-I--I- rewrite of the original 
Fortran Minuit package [42] . 


C. Estimating uncertainties of PDFs 


In section III HI we described how we obtain our best 
estimate (the central value) of the nCTEQ nuclear PDF s 
as the minimum of the x^ function defined in Eq. (2.121. 
Now we want to probe the vicinity of this minimum to 
be able to estimate uncertainties on our prediction. This 


is done using the Hessian method [331 Sll j which will be 
briefly described in the following. We follow the notation 
of Ref. [33] and refer the reader to this publication for 
more details on the Hessian formalism. 


1. Determination of the Hessian matrix 


The basic assumption of the Hessian method is that 
near its minimum the x^“function can be approximated 
by a quadratic form of the fitting parameters {ai}. 
Therefore, it can be written as 


X^ = xl^Y^Hijy,yj, (2.16) 


where yi = ai — are the parameter shifts from the 
minimum given by the a° parameters, Xq = X^({Qi}) is 
the value of the x^-function in the minimum, and Hij is 
the Hessian matrix defined as: 


H^J = 


1 

2 


dyidyj 



(2.17) 


Since the Hessian Hij is a symmetric nxn matrix (where 
n is the number of free parameters ai) it has n orthogo¬ 
nal eigenvectors forming a basis in the {?/i}“Space. The 
characteristic equation can be written as: 

(2.18) 

3 


(k) 

The eigenvectors ^ that we use can be normalized so 
that: 


(2.19) 


For our later convenience we also introduce eigenvectors 
normalized to the corresponding eigenvalues: 


y(fe) 

I 


1 


yik) 

I 


( 2 . 20 ) 


The eigenvectors can be used to disentangle the original 
PDF parameters and define a new basis z = {zi} where 
the Hessian is diagonal^ 


Y. y^ = 

‘i-J 

( Ai 

j' 0 

= z . 

V 0 


Zj = z^.D'^.H.D.z 


0 

A 2 


0 \ 


0 

A„ / 


( 2 . 21 ) 


^ In the basis defined using the rescaled eigenvectors the 

Hessian is represented by a unit matrix. 
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The new coordinates are defined using a matrix D as 

z = i?-V, (2.22) 

where D is & matrix composed of eigenvectors: 





. \ 




. ^ 


VU« 


. ) 


(2.23) 


Note that because the Hessian is symmetric D ^ 

(i) 

Using the index notation such as Dij = we can write 
the relation between the original fitting parameters and 
the new parameters as: 


Vi = 




E 




(2.24) 


where we introduced a new basis Zi which corresponds to 

*■ (k) 

the rescaled eigenvectors 1^ . The inverse transforma¬ 

tion is given by: 


j 

3 3 


(2.25) 


In the new coordinates, ~Xo ^ particularly 

simple form: 


= . (2.26) 
i i 

Using the Hessian method to analyze the vicinity of 
the minimum of the y^-function seems straightforward 
in theory but in practice when applied to a global PDF 
analysis, one encounters a few problems worth pointing 
out. As was already mentioned in the discussion of free 
proton PDFs [13] and as is the case in our analysis, the 
eigenvalues of the Hessian span several orders of magni¬ 
tude. In order to correctly identify all eigenvalues, the 
precision with which the Hessian matrix is determined 
needs to be kept under control. 

In practice, the Hessian matrix is calculated using fi¬ 
nite differences to determine the second derivatives. A 
careful choice of the step in the finite difference defini¬ 
tion of the second derivatives is crucial. If the step is 
too large, one probes too large a neighborhood of the 
minimum where the y^-function cannot be described by 
a quadratic approximation anymore. If the step is too 
small, numerical noise in the y^-function prevents a reli¬ 
able determination of the second derivatives. Moreover, 
the step size has to be different for each of the parameters 
as the y^-function depends differently on each of them. 


The relative step sizes Ayi to each of the parameters are 
set as 


Ayi = 



(2.27) 


where Ay^ = x^ ~ Xo defines the small neighborhood 
from which the derivatives of the y^-function are calcu¬ 
lated. 

It turns out that the numerical noise in the y^-function 
is larger than expected for the case of a global PDF analy¬ 
sis. Contrary to what one would expect, the y^-function 
is not smooth which influences the determination of the 
second derivatives for all step sizes. It all comes down to 
the fact that one evaluation of the y^-function requires 
several hundred evaluations of different next-to-leading 
order theory calculations which, in their numerical im¬ 
plementations, are not smooth functions of the fit pa¬ 
rameters. 

To reduce the influence of the noise on the derivatives 
of the y^-function, the standard definition of the deriva¬ 
tive using the central differences 


df _ U, - /_! 
dx 2h 


(2.28) 


in which fk = f{xo + kh), is replaced by noise reducing 
derivatives (see [45] ). The central differences approach to 
derivatives is based on interpolating the y^-function by a 
polynomial which coincides with the y^-function in sev¬ 
eral chosen points e.g. a quadratic polynomial interpo¬ 
lating t he y^ -function in 3 points leads to the derivative 
in Eq. (2.28). If the y^-function suffers from numerical 
noise, the interpolated polynomial suffers as much if not 
more. 

We adopt a different approach and instead of interpo¬ 
lating N points by a polynomial of the order A — 1, we 
allow a polynomial to assume different values in these N 
points and approximate the y^-function by the method 
of least squares. This approach assumes that the order of 
the polynomial M has to be strictly less than N—1 where 
N is the number of points. If we use a quadratic poly¬ 
nomial to fit 7 symmetrically chosen, equidistant points 
of the y^-function, we obtain the following prescriptions 
for the 7-point low-noise derivative 


df _ fi — f-i + 2(/2 — f- 2 ) + 3(/3 — /_ 3 ) 

dx~ 28h ■ i ■ ) 


Using these derivatives instead of the standard derivative 
from Eq. (2.28) and extending this approach to the sec¬ 
ond derivatives allows us to determine the Hessian with 
sufficient precision and to eliminate the influence of the 
numerical noise. 


2. Error PDFs 

To translate the uncertainties contained in the data 
to the underlying PDF parameters, we use the fact that 
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the x^-function in the diagonalized Hessian approxima¬ 
tion is a simple function of the parameters Zk- Varying 
data within their errors corresponds to a change in 
(denoted by Ay^) which can then in turn be interpreted 
as a shift in the parameters Zk 

Zk = ±VAy2, 

[a^ (2.30) 

A specific change in y^ can be obtained by varying the 
parameters using n independent directions in the param¬ 
eter spacej^ In the Zk space all directions are equivalent 
so we can choose the n independent directions to coincide 
with the directions where one single parameter is varied. 
A change in one direction along one single parameter Zk 
leads to a simultaneous change in all original parameters 

Qi 


2/. = Aa, = ±y^v/"). (2.31) 

The parameter shifts along the direction of the Zk param¬ 
eter are used to generate 2n error PDFs for a specified 
Ay2 


fk^f 





for k = 1, 2,..., n. 


(2.32) 


The error PDFs can be used to determine the PDF un¬ 
certainty of any observable X which depends on PDFs. 
This uncertainty, which we denote as AA, can be deter¬ 
mined in different ways and in this work we define it by 
adding errors in quadrature 

AV = iXifk)-Xifk))" . (2.33) 


The PDF uncertainty AX clearly depends on the exact 
value chosen for Ay^. In an ideal case, an increase of y^ 
corresponding to one standard deviation from the central 
value is Ay^ = 1. However, in our fit we combine results 
from different experiments which are not necessarily un¬ 
correlated or compatible, so the standard argument does 
not apply and Ay^ may be different from one. To esti¬ 
mate what is the appropriate value for the Ay^ (often 
referred to as the tolerance) we use a criterion similar to 
the one advocated in miisiiiz], which results in the 
value Ay^ = 35. Additionally, since the value of our 
tolerance is far from 1, the quadratic approximation of 


If one allows only positive changes of parameters, there are 2n 
directions. 


the Hessian method becomes less precise. We account 
for it by introducing an additional procedure of rescaling 
of the Hessian matrix. Both the rescaling procedure and 
the criterion for choosing the Ay^ tolerance are described 
in detail in Appendix [A| 


III. EXPERIMENTAL DATA 


F^/F° : 
Observable 

Experiment 

ID 

Ref. 

# data 

^ data 

after cuts 


D 

NMC-97 

5160 

EH 

292 

201 

247.73 

He/D 

Hermes 

5156 

EH 

182 

17 

13.45 


NMC-95,re 

5124 

Efll 

18 

12 

9.78 


SLAC-E139 

5141 

IsTl 

18 

3 

1.42 

Li/D 

NMC-95 

5115 

EH 

24 

11 

6.10 

Be/D 

SLAC-E139 

5138 

EH 

17 

3 

1.37 

C/D 

FNAL-E665-95 

5125 

1531 

11 

3 

1.44 


SLAC-E139 

5139 

ED 

7 

2 

1.36 


EMC-88 

5107 

El 

9 

9 

7.41 


EMC-90 

5110 

ED 

9 

0 

0.00 


NMC-95 

5113 

El 

24 

12 

8.40 


NMC-95,re 

5114 

Isol 

18 

12 

13.29 

N/D 

Hermes 

5157 

El 

175 

19 

9.92 


BCDMS-85 

5103 

l56l 

9 

9 

4.65 

Al/D 

SLAC-E049 

5134 

m 

18 

0 

0.00 


SLAC-E139 

5136 

El 

17 

3 

1.14 

Ca/D 

NMC-95,re 

5121 

l50l 

18 

12 

11.54 


FNAL-E665-95 

5126 

l53l 

11 

3 

0.94 


SLAC-E139 

5140 

El 

7 

2 

1.63 


EMC-90 

5109 

ED 

9 

0 

0.00 

Fe/D 

SLAC-E049 

5131 

El 

14 

2 

0.78 


SLAC-E139 

5132 

El 

23 

6 

7.76 


SLAC-E140 

5133 

l59l 

10 

0 

0.00 


BCDMS-87 

5101 

El 

10 

10 

5.77 


BCDMS-85 

5102 

El 

6 

6 

2.56 

Cu/D 

EMC-93 

5104 

El 

10 

9 

4.71 


EMC-93 (chariot) 

5105 

El 

9 

9 

4.88 


EMC-88 

5106 

El 

9 

9 

3.39 

Kr/D 

Hermes 

5158 

El 

167 

12 

9.79 

Ag/D 

SLAC-E139 

5135 

El 

7 

2 

1.60 

Sn/D 

EMC-88 

5108 

El 

8 

8 

17.20 

Xe/D 

FNAL-E665-92 

5127 

El 

10 

2 

0.72 

Au/D 

SLAC-E139 

5137 

El 

18 

3 

1.74 

Pb/D 

FNAL-E665-95 

5129 

El 

11 

3 

1.20 

Total: 




1205 

414 

403.70 


Table I: The DIS F^/F.^ data sets used in the nCTEQlS 
fit. The table details values of y^ for each experiment, 
the specific nuclear targets, references, and the number 
of data points with and without kinematic cuts. 


In the current analysis we use deep inelastic scattering 
data (DIS), Drell-Yan lepton pair production data (DY) 
and inclusive pion production data from RHIC (for nuclei 
with A > 2). The details of particular experiments such 
as the number of data points, measured observables, etc. 
are summarized in Tables Ul lIVI 

The reason to include data from different processes is 
























F*/F*' : 
Observable 

Experiment 

ID 

Ref. 

^ data 

^ data 

after cuts 


c/Li 

NMC-95,re 

5123 

I50l 

25 

7 

5.56 

Ca/Li 

NMC-95,re 

5122 

l50l 

25 

7 

1.11 

Be/C 

NMC-96 

5112 

m 

15 

14 

4.08 

Al/C 

NMC-96 

5111 

111] 

15 

14 

5.39 

Ca/C 

NMC-95,re 

5120 

[50] 

25 

7 

4.32 


NMC-96 

5119 

m] 

15 

14 

5.43 

Fe/C 

NMC-96 

5143 

l63l 

15 

14 

9.78 

Sn/C 

NMC-96 

5159 

l6ll 

146 

111 

64.44 

Pb/C 

NMC-96 

5116 

m] 

15 

14 

7.74 

Total: 




296 

202 

107.85 


Table II: The DIS F.^ jF^ data sets used in the 
nCTEQlS fit. We list the same details for each data set 
as in Tab. 



pA j pA' 
^DV/'^DV • 
Observable 

Experiment 

ID 

Ref. 

^ data 

^ data 

after cuts 


C/H2 

FNAL-E772-90 

5203 

I65l 

9 

9 

7.92 

Ca/H2 

FNAL-E772-90 

5204 

[65l 

9 

9 

2.73 

Fe/H2 

FNAL-E772-90 

5205 

l65l 

9 

9 

3.17 

W/H2 

FNAL-E772-90 

5206 

l65l 

9 

9 

7.28 

Fe/Be 

FNAL-E886-99 

5201 

[66] 

28 

28 

23.09 

W/Be 

FNAL-E886-99 

5202 

[66] 

28 

28 

23.62 

Total: 




92 

92 

67.81 


Table III: The Drell-Yan process data sets used in the 
nCTEQlS fit. We list the same details for each data set 
as in Tab. 


^dAu/^pp ■ 

Observable 

Experiment 

ID 

Ref. 

^ data 

^ data 

after cuts 


dAu/pp 

PHENIX 

PHENIX 

EZ] 

21 

20 

6.63 


STAR-2010 

STAR 

[H 

13 

12 

1.41 

Total: 




34 

32 

8.04 


Table IV: The pion production data sets used in the 
nCTEQlS fit. We list the same details for each data set 
as in Tab. 


Figure 1: Kinematic reach of DIS and DY data used in 
the presented nCTEQ fits. The dashed lines represent the 
kinematic cuts employed in this analysis (Q > 2 GeV, 
W > 3.5 GeV). Only the data points lying above both 
of these lines are included in the fits. 



X 

Figure 2: Approximate x-range for the pion data with 
the Binnewies-Kniehl-Kramer fragmentation function. 


that each process helps constrain different combinations 
of parton distributions. The bulk of our data are from 
DIS which help pin down the valence and sea distribu¬ 
tions, however they are not very sensitive to different 
quark flavors and gluons. The DY data can be used to 
differentiate between u and d quark flavors, and the in¬ 
clusive pion data have a potential to better constrain the 
gluon distribution]^ 


® Note that the inclusive pion production observable is different in 
the sense that it has an additional dependence on a fragmentation 
function. 


We introduce kinematic cuts on the included data 
which limit possible effects of higher twist contributions 
and target mass corrections and at the same time are 
compatible with the kinematic cuts used in the underly¬ 
ing free proton analysis. The cuts used in this analysis 
are: 

• DIS: Q > 2 GeV and W > 3.5 GeV, 

• DY: 2 < M < 300 GeV, 

(where M is the invariant mass of the produced 
lepton pair) 

• 7r° production: py > 1.7 GeV. 
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After the cuts are applied, 740 data points remain, 
including 616 DIS, 92 DY and 32 pion production data 
points. 

Note that the overall number of data points we use is 
considerably smaller compared to the number of data fit¬ 
ted by other groups (e.g., EPS [12] has 929 data points). 
One reason is that the other analyses employ less strin¬ 
gent kinematic cuts on Q^: 

• EPS [I2|: Q > 1.3 GeV, 

• HKN mj: Q > 1 GeV, 

• DSSZ [U: Q > 1 GeV. 


In addition, none of the analyses mentioned above em¬ 
ploy a cut on W. Whereas the looser cuts allow one to 
use more data in the fit, there are possible disadvantages 
connected to this choice. In particular, if one adopts loose 
cuts, one runs into the danger that the contributions from 
the target mass effects or higher twist effects can get en¬ 
hanced. Especially the latter effects may be more im¬ 
portant in the nuclear case due to the higher density of 
spectator partons in the nucleus usmsi and so their ef¬ 
fect can be easily underestimated. However, the effect of 
higher twist and target mass corrections have been shown 
to be weakened in ratios of observables [691 ED]. 

The kinematic reach of the DIS and DY data sets used 
in our fit is summarized in Fig.j^ where individual exper¬ 
imental points are shown in the (x, Q^) plane. Note that 
the two dashed lines indicate the kinematic cuts; points 
lying below these lines are excluded from our analysis. 

In Fig. we estimate the kinematic impact of the 
pion data by plotting the cross section for inclusive pion 
prod ucti on before convoluting it with gold PDFs, see 
Eq. (2.4). Fig. shows the normalized cross-section as 
a function of the Bjorken-x of a parton inside a nucleon 
of a gold atom. This is only an estimate which uses the 
leading order (LO) prediction and it also depends on the 
fragmentation function (FF) that is used. Nevertheless, 
it is useful and allows us to see that the x-values probed 
by the pion data depend quite substantially on the pr- In 
particular, for higher px, higher x values are probed, e.g. 
for pt ~ 15 GeV we are mostly sensitive to x G (0.2,0.3), 
whereas for lower px the probed x values are more dif¬ 
fused, e.g. for Pt ~ 2 GeV x G (0.01,0.04). 

One should mention that there are still experimental 
data that could have been included in our analysis but we 
have decided for different reasons to exclude them from 
the current work. We comment briefly on the two most 
important examples. 

First, there are neutrino DIS data from CDHSW in], 
GHORUS [72], and in particular from the NuTeV collab¬ 
oration m- Since they include a considerable number 
of data points and probe more flavor combinations than 
the charged lepton data, they can be used to differen¬ 
tiate individual flavors. However, tensions between the 
inclusive charged current i^A DIS data from NuTeV and 
the neutral current i^A data found in [nnniiii] indi¬ 
cate that some additional effort is required to understand 


how these discrepancies can be resolved so that all data 
could be used in one fit simultaneously. Since these dis¬ 
crepancies appear only if one takes into account the full 
information contained in the correlated error matrix, ne¬ 
glecting these correlations makes it possible to combine 
jyA and ^^A DIS in one fit [TTl l22l [74] . We plan to revisit 
the neutrino data in a future publication but decided not 
to include them in our present PDF release. 

Another important set of data which could be included 
are the already available LHG data. In particular the 
cleanest probe of nuclear effects at the LHC comes from 
the vector boson, W^, Z, production [7M7H]- Results 
on asymmetries in pPb collisions |78j in particular have 
a potential to provide valuable input for nuclear PDF 
analyses. These data are not included in the current 
release as we first want to provide a baseline analysis 
without any LHC data. 


IV. RESULTS 


A key result of the current nCTEQlS fit compared to the 
previous nCTEQ releases [HED] is the inclusion of PDF 


uncertainties using the Hessian method, cf. Sec. H C 


The second significant addition is the inclusion of a 
new type of experimental data, namely the pion pro¬ 
duction data from the PHENIX and STAR collabora¬ 
tions. Since these data have the potential to provide 
information on the gluon distribution (which otherwise 
is weakly constrained) it is important to precisely esti¬ 
mate their impact on the resulting PDFs. For this pur¬ 
pose the nCTEQlB fit will be compared with a reference fit 
nCTEQ15-np which is identical except it does not include 
the pion data. 

The full set of data we consider is listed in 
Tables IIHI3 Note that we have included QED radiative 
corrections for the DIS FNAL-E665-95 (Pb/D, Ca/D, 
C/D) data sets and this significantly improves the de¬ 
scription of these dataj^ In the following, we discuss 
the results of this nCTEQlB analysis and compare it with 
other available sets of nuclear PDFs. 


A. The nCTEQlS Fit 


1. PDF Parameterization 


The PDFs in our fit are parameterized at the input 
scale Qo = 1-3 GeV according to Eqs. (2.5) and (2.6). 
This provides considerable flexibility as each of the seven 
flavor combinations can have ~10 free parameters to de- 


For example, the for the FNAL-E665-95 Pb/D data (ID 5129) 
is reduced from 5.91 to 1.20 (for 3 data points) when the QED 
radiative corrections are included. 
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Par. 

Value 

Par. 

Value 

Par. 

Value 

Par. 

Value 

Par. 

Value 

Par. 

Value 

M‘> 

(0.382) 


(0.327) 


(0.136) 

jy^d+u 

(0.129) 

M‘>+^ 

(0.026) 

Md/u 

(0.000) 

^0,0 

(0.000) 


- 

- 

- 

- 

- 

^s+s 

^0,0 

(0.500) 

- 

- 

^1,0 

(0.523) 

^1,0 

(0.630) 

^1,0 

(0.513) 

^1,0 

(-0.324) 

^s+s 

^1,0 

(-0.324) 

diu 

^1,0 

(10.075) 

^2,0 

(3.034) 

^2,0 

(2.934) 

^2,0 

(4.211) 

d+u 

^2,0 

(8.116) 

^2,0 

(8.116) 

dju 

^2,0 

(4.957) 

^3,0 

(4.394) 

^3,0 

(-2.369) 

^3,0 

(-2.375) 

d+u 

^3,0 

(0.413) 

^s+s 

^3,0 

(0.413) 

d/u 

(15.167) 

^4,0 

(2.359) 

^4,0 

(1.266) 

^dv 

^4,0 

(0.965) 

d+u 

^4,0 

(4.754) 

^4,0 

(4.754) 

d/u 

^4,0 

(17.000) 

^5,0 

(-3.000) 

^5,0 

(1.718) 

^5,0 

(3.000) 

d+u 

^5,0 

(0.614) 

^s+s 

^5,0 

(0.614) 

d/u 

^5,0 

(9.948) 

Par. 

Value 

Par. 

Value 

Par. 

Value 

Par. 

Value 

Par. 

Value 

Par. 

Value 

^0,1 

(-0.256) 

- 

- 

- 

- 

- 

- 

^s+s 

^0,1 

(0.167) 

- 

- 

^ 1,1 

- 0.001 

^1,1 

- 2.729 

d^ 

^1,1 

0.272 

d+u 

^1,1 

0.411 

^S + S 
^1,1 

(0.411) 

d/u 

^1,1 

(0.000) 

^2,1 

(0.000) 

14,, 

^2,1 

- 0.162 

dv 

^2,1 

- 0.198 

d+u 

^2,1 

(0.415) 

^S + S 
^2,1 

(0.415) 

d/u 

^2,1 

(0.000) 

^3,1 

(0.383) 

^3,1 

(0.018) 

^3,1 

(0.085) 

d+u 

^3,1 

(-0.759) 

^S + S 

^3,1 

(0.000) 

d/u 

^3,1 

(0.000) 

^4,1 

0.055 

^4,1 

12.176 

^d^ 

M,! 

(3.874) 

d+u 

^4,1 

(-0.203) 

^S + S 

^4.1 

(0.000) 

d/u 

^4,1 

(0.000) 

^5,1 

0.002 

^5,1 

- 1.141 

f.dv 

^5,1 

- 0.072 

d+u 

^5,1 

- 0.087 

^S+S 

^5,1 

(0.000) 

d/u 

^5,1 

(0.000) 

Par. 

Value 

Par. 

Value 

Par. 

Value 

Par. 

Value 

Par. 

Value 

Par. 

Value 

^9 

^0,2 

- 0.037 

- 

- 

- 

- 

- 

- 

^s+s 

^0,2 

(0.104) 

- 

- 

^ 1,2 

- 1.337 

^ 1,2 

(0.006) 

^ 1,2 

(0.466) 

d+u 

^1,2 

(0.172) 

^s+s 

^ 1,2 

(0.172) 

d/u 

‘- 1,2 

(0.000) 

^ 2,2 

(0.000) 

^ 2,2 

(0.524) 

^ 2,2 

(0.440) 

d+u 

^2,2 

(0.290) 

^s + s 
^2,2 

(0.290) 

d/u 

<- 2,2 

(0.000) 

f.9 

^ 3,2 

(0.520) 

^ 3,2 

(0.073) 

^ 3,2 

(0.107) 

d+u 

^ 3,2 

(0.298) 

^s + s 
^ 3,2 

(0.000) 

d/u 

^ 3,2 

(0.000) 

^ 4,2 

- 0.514 

^ 4,2 

(0.038) 

oia 

(-0.018) 

d+u 

^ 4,2 

(0.888) 

^ 4,2 

(0.000) 

d/u 

<- 4,2 

(0.000) 

f.9 

^ 5.2 

- 1.417 

^ 5.2 

(0.615) 

0^2 

(-0.236) 

d+u 

^ 5.2 

(1.353) 

^ 5.2 

(0.000) 

d/u 

^ 5.2 

(0.000) 


Table V: Values of the parameters of the nCTEQ15 fit at the initial scale Qo = 1-3 GeV. Values in bold represent the 
free parameters and values in parentheses are fixed in the fit. The not listed normalization parameters are 
determined by the momentum and number sum rules as discussed in the text. For completeness, we provide the full 
set of the free proton parameters Ckfi (first set of rows). The M* parameters (first row) show the (fixed) momentum 

fraction carried by different flavors in the case of a free proton. 


scribe the x and A dependence!^ However, the available 
experimental data are not sufficient to constrain such a 
flexible parameterization. Therefore, we limit our actual 
fit to 16 parameters; specifically, we include 7 gluon, 4 
M-valence, 3 d-valence and 2 d + u free parameters. The 
details of the fit are summarized in Table IVl which shows 
the best fit values of the free parameters, as well as the 
values of the fixed parameters. 

For the pion data, we allow for the normalization to 
vary and we obtain 1.031 for the PHENIX data m and 
0.962 for the STAR data [5H]|^ Our obtained normal¬ 
ization shifts of ~4% lie well within the experimental 
normalization uncertainty!^ 

Our parameterization smoothly interpolates between 


For each of the 5 flavor combinations -|- d, s = s} 

of Eq. ( |2.5| l we have 10 parameters {cj, i,c^. 2 } for k = {1...5} 
in addition to the normalization parameters cq that are partly 
fixed by the number and momentum sum rules. For {d/u}, we 
have 8 parameters at our disposal. 

Note that the data normalization parameters do not enter the 
Hessian analysis of uncertainties. 

See Table 1 and Fig. 2 in Ref. [57] for PHENIX, and Fig. 25 in 
Ref. 1681 and Table 5 in Ref. 1791 for STAR. 


different nuclei as a function of the nuclear mass number 
A] the number of protons Z and neutrons {A — Z) en¬ 
ters only through the isospin composition of a nucleus, 
cf. Eq. (2.1). Fig. !^ shows the A dependence of the fit¬ 
ting parameters normalized by the corresponding values 
of the free proton baseline parameters Ck o- (Note, some 
of these parameters are fixed, cf. Table [^ ) Many of the 
parameters change rapidly in the region of light nuclei 
A < 25 and are relatively stable for heavy nuclei A > 50. 
Also, we observe that the parameters responsible for the 
small X behavior {ci}, typically exhibit a strong A depen¬ 
dence, whereas the large x parameters { 02 } are compara¬ 
bly insensitive to the type of nucleus. In particular, the 
biggest effect occurs for the gluon where the cf param¬ 
eter describing the low x gluon PDF and c| parameter 
(responsible for mid-cc) are changing linearly throughout 
the whole range of A. 


2. of Jd 

We now examine the overall statistical quality of the 
fit as measured by the For the nCTEQlS fit we obtain 
a total of 587.4 with 740 data points (after kinematic 
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(a) Gluon 



(c) u-valence 



(b) d-\- u 



Figure 3: ^-dependence of the fit parameters as given in Eq. (2.6). Specifically, we plot 
Ck{A) = Ck,o + Cfcu(l — for each flavor normalized to the corresponding free proton para mete r Ck,o- The 

superscripts {1, 2, ...} in the legend correspond to the parameters {ci, C 2 , ... } in Eq. (2.5). 


cuts). With 18 free parameters (including 2 data normal¬ 
ization parameters) this leads to a x^/do/ = 0.81 which 
indicates a good fit. Furthermore, this /dof is not too 
small which could indicate deficiencies of the fit such as 
over-fitting. 


To better evaluate the fit quality, in Fig. we plot 
the /dof for the individual experiments and check that 
the majority of experiments has a {x^/dof) ~ 1. While 
most experiments satisfy this “goodness of fit” criterion, 
there is one experiment that stands out as having a poor 
fit: the DIS EMC-88 data for Sn/D (ID 5108). Several 
previous global analyses have also found it challenging to 
accommodate the Sn/D data [TTIIH] . 


In Fig. 4b we show again the x^/do/, but this time the 
experiments are grouped by nuclear target and are sorted 
by increasing nuclear mass number A. This allows us to 
see that there are no systematic effects associated with 
our choice of the A parameterization. With the noted 
exception of Sn/D, all other nuclear targets from helium 
up to lead are described very well with a ^/dof ~ 1. 


3. Error PDF reliability. 

Before we examine the actual nCTEQlS predictions, we 
first investigate the quality of the Hessian error analysis. 
This will allow us to judge the reliability of our error 
estimates and, in turn, the quality of our predictions!^ 
There are two factors that need to be assessed: 

(i) the quality of the quadratic approximation, 

(ii) how well the Hessian approximation describes the 
actual x^ function in a region around the minimum 
given by our tolerance criterion, = 35. 

To estimate these factors, we plot the x^ function relative 
to its value at the minimum (Ay^ = X^~Xo) along the 16 


Note that by construction, the Hessian method can only probe 
the local minimum connected to the “best fit” (central predic¬ 
tion), and is not sensitive to a landscape with multiple minima. 
Unfortunately, in case of nPDFs fits multiple minima are possible 
as there is not sufficient data to fully constrain the nPDFs. 
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(a) Value of x^/dof for the individual experiments which are identified by the IDs that are listed in 
Tables m-flVl The 51xx IDs correspond to the DIS experiments, the 52xx IDs are the DY data, and the pion 
data are labeled by the collaboration name. The experiments are sorted left-to-right: {DIS, DY, tt®} and 

sub-sorted by the nuclear mass number A. 



(b) Value of x^/dof per nuclear target used in the nCTEQlS fit sorted left-to-right by the nuclear mass number A. 


Figure 4: Value of x^/dof for (a) individual experiments and (b) per nuclear target used in the nCTEQlS fit. The 
numbers on top of the bars represent the number of data points (after kinematic cuts). 


error directions in the eigenvector space (see Fig.[^. For 
comparison, we also display the Hessian approximation 
given by the quadratic form = zf. The plots are or¬ 
dered according to the decreasing values of the eigenval¬ 
ues corresponding to the Zi directions; the largest eigen¬ 
value is of order 10®, and the smallest of order 10. For the 
largest few eigenvalues of Fig. the quadratic approxi¬ 
mation works extremely well; however, for the smaller 
eigenvalues {e. 3 ., #10, #14} it can deviate from the 
function. Nevertheless, in all the cases we are able to 
obtain a good description of the actual y® function for 
Zi ~ [— 6 , 6 ] which corresponds to our tolerance criterion 
yAy® = -#35 ~ 6 . This analysis verifies that the error 
PDFs defined using the modified Hessian formalism will 
closely reflect the actual y® function determined by the 
experimental data, and will not be severely affected by 
the imperfections of the quadratic approximation that oc¬ 


curs for directions corresponding to lower eigenvalues 


4- nPDFs vs. nuclear A 

We now examine the results of the nCTEQlB fit starting 
with the H-dependence of the various nPDF flavors. In 
Fig. [^we display the central fit predictions for a range of 
nuclear A values from A = 1 (proton) to A = 208 (lead). 
When examining the A-dependence we observe that as 
we move to larger A the gluon and sea-distributions 
{g^u^d^s} decrease at small x values. This trend is also 


In the modified Hessian approach that we use, the discrepancies 
at Ax^ = 35 originate mostly from the non-symmetric behavior 


of the x^ function; see Sec. 


lie 


and 


Appendix [a| 


for details. 
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Figure 5: function relative to its value at the 

minimum, Ax^ = ~ Xoi plotted along the 16 error 

directions in the eigenvector space, zf. We display the 
true x^ function (solid lines) and the quadratic 
approximation given by Hessian method = zf 
(dashed lines). The eigenvector directions are ordered 
from the largest to the smallest eigenvalue. 


present for the {u, d} PDFs. On the other hand, the A- 
dependence of {uy,dv\ distributions is reduced relative 
to the other flavor cornponents. 

Finally, Figs. 0 and§ show our nPDFs for a 

lead nucleus together with the nuclear correction factors 
at the input scale Q = Qo = \.2> GeV and at Q = 10 GeV 
to show the evolution effects when the PDFs are probed 
at a typical hard scale. We have chosen to present results 
for the rather heavy lead nucleus because of its relevance 
for the heavy ion program at the LHG. In all cases, we 
display the uncertainty band arising from the error PDF 
sets based upon our eigenvectors and the tolerance crite¬ 
rion. It should be noted that the uncertainty bands for 
^ ^ 10”^ and X > 0.7 are not directly constrained by 
data but only by the momentum and number sum rules. 
The uncertainty bands are the result of extrapolating the 
functional form of our parametrization into these uncon¬ 
strained regions. 

Some comments are in order: 

• As can be seen from Fig. (a), our input gluon is 
strongly suppressed/shadowed with respect to the 
free proton in the x < 0.04 region. In fact, it has a 
valence-like structure (see Fig. (b)) which van¬ 
ishes at small x. Gonsequently, the steep small 
X rise of the gluon distribution at Q = 10 GeV 
(see Fig. is entirely due to the QCD evolution. 





Figure 6: nCTEQlS bound proton PDFs at the scale 
Q = 10 GeV for a range of nuclei from the free proton 
{A = 1) to lead {A = 208). 


However, we should note that there is no data con- 
strints below x ~ 0.01 and the gluon uncertainty 
in this region is underestimated. In addition, our 
gluon has an anti-shadowing peak around x ~ 0.1 
and then exhibits suppression in the EMC region 
X ^ 0.5. However, the large x gluon features wide 
uncertainty band reflecting the fact that there are 
no data constraints. 

• In our analysis we determine the u + d combination 
and assume that there is no nuclear modification 
to the d/u combination (see Sec. [lT| and Table 0- 
As a result the u and d PDFs are very similar, the 
small difference between the two comes from the 
underlying free proton PDFs. 

• In this analysis we do not fit the strange distribu¬ 
tion but relate it to the light quarks sea distribu¬ 
tion, see Eq. (2.7). As a result the strange quark 


distribution is very similar to the u and d distribu¬ 
tions. 
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Figure 7: Results of the nCTEQlS fit. On the left we show nuclear modification factors dehned as ratios of proton 
PDFs bound in lead to the corresponding free proton PDFs, and on the right we show the actual bound proton 
PDFs for lead. In both cases the scale is equal to Q = 1.3 GeV. 


• Contrary to the other existing nPDFs where the 
nuclear correction factors for the valence distribu¬ 
tions are assumed to be the same, we treat and 
dy as independent. This leads to an interesting 
feature of our result where Uy is suppressed and 
dy is enhanced in the EMC region. This behavior 
is not entirely unexpected, there are nuclear mod¬ 
els predicting a flavor dependence for the EMC ef¬ 
fect [Tomniisi]- 

• The above difference for the nuclear correction in 

Uy and dy appears at the level of the bound proton 
PDFs. When we construct a physical combination 
representing the full nuclear PDF, -I- 

such as lead in Fig. 9 the combination 
yields net corrections for Uy anoa^ which are close 
to each other and similar to those in the literature. 
We will discuss this in more detail in Sec. IIV El 

Once more data are included, e.g. from the LHC, 


neutrino DIS experiments and a future eA collider, 
it should be possible to relax some of the assump¬ 
tions. 

In the following section, we will investigate the impact 
of these nPDFs and the corresponding uncertainty bands 
on the physical observables. 


B. Comparison with data 

While the y^/dof is one measure of the quality of the 
fit this alone obviously does not capture all the relevant 
characteristics. To investigate the nCTEQlB result in more 
detail we compare it to the most important and con¬ 
straining data sets and consider strengths and limitations 
of both the ht and the available data sets. 
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Figure 8: Results of the nCTEQlB fit. On the left we show nuclear modification factors dehned as ratios of proton 
PDFs bound in lead to the corresponding free proton PDFs, and on the right we show the actual bound proton 
PDFs for lead. In both cases the scale is equal to Q = 10 GeV. 


1. DIS data sets 


The data from the deep-inelastic scattering experi¬ 
ments are by far the most numerous and provide the 
dominant contribution to the total These experi¬ 
ments are performed on a variety of nuclei which allow us 
to constrain the A dependence of our parameters. Most 
of the data are extracted as a ratio of structure func¬ 
tions R = I for two different targets A\ and A 2 . 
Note that in the present study we do not fit data from 
the very high x region x > 0.7 since they do not pass 
our kinematic cuts. As already mentioned the high x re¬ 
gion is theoretically challenging due to a host of effects 
(higher twist, target mass corrections, large x resumma¬ 
tion, deuteron wave-function, nuclear off-shell effects). 
Some of these effects in the large x and low area have 
been investigated extensively in the proton case by the 
CTEQ-CJ collaboration [5HIH3]. The nuclear case is even 
more challenging due to enhanced higher twist and Fermi 


motion effects which lead to a steep rise of the structure 
function ratios in the limit a; —)■ 1. For these reasons we 
avoid fitting the high x region for the time being. 

The comparison of our fit to the DIS F 2 ratio data is 
shown in Figs. and 11 as a function of x. Note, in 
these figures the data for different are combined into 
a single plot as the scaling violations (discussed later) 
occur on a logarithmic scale and largely cancel out in the 
ratios. 

shows the ratio F^{x, Q^)lF^{x, Q^) for a va¬ 


Fig. 


10 


riety of experiments. The overall agreement of the fit 
with the data is excellent for a majority of the nuclei. 
The discrepancy which can be seen for the EMC data 
taken on tin (Sn/ D) is t he same discrepancy we have 
pointed out in Sec. IV A 2 when we investigated the of 
the individual experiments. As already mentioned, this 
problem has been also encountered in previous analyses 
[in HU and we are unable to reconcile it with our ht. 

Similarly, Fig. shows the structure function ratio 
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Figure 9: We show nuclear modification factors defined as ratios of lead PDFs compared to a lead PDF constructed 
from free-proton PDFs. The PDFs are constructed using for ^°^Pb and the free proton, 

with a scale of Q = 10 GeV. 


F^{x, (x, in comparison to NMC data for a 

variety of nuclear targets. These high-statistics data are 
also well described by the results of the nCTEQlS fit. 

The NMC data taken on tin and carbon (i? = 
F^"^ I f!?) c over a wider range in and we display these 

in Fig. as a function of binned in x. As is well 
know, the logarithmic scaling violations of the struc¬ 
ture functions provide constraints on the low x gluon dis¬ 
tribution. Of course, compared to the very precise HERA 
data on the proton F 2 structure function which extends 
over a very wide range of values the NMC data have a 
much smaller lever arm. As a consequence the NMC 
data provide relatively weaker constraints on the nuclear 
gluon PDF in the x range of (0.05,0.1). We will discuss 
data constraints on gluon in more detail in Sec. |IVD[ 


In Fig. 13 we plot the nuclear correction R = Ff^/F^ 
for iron vs. x for two values and compare the results 
with experimental data and with results from different 
nPDF groups. Comparing these two figures, we again 
see that there is a rather weak Q^-dependence of the 
structure function ratio between = 5 GeV^ and = 
20 GeV^. As discussed above due to our strict kinematic 
cuts we do not extend our predictions to the high x region 
{x > 0.7). 


Taking into account both the nPDF uncertainty (rep¬ 
resented by the error bands) and the experimental error 
bars, the data are generally compatible with the nCTEQlS 
fit. In addition to comparing with data, we compare our 
predictions with those of HKN [2] and EPS [2] and find 
a good agreement within the errors of our analysis. 


2. Drell-Yati data sets 


We now turn to the Drell-Yan muon pair produc¬ 
tion process p + A^^'^+^~+X. In Fig. 14 (a). 


we display the differential cross section ratio, R = 
{da^/dx 2 dM)/{da^/dx 2 dM), measured by the Fer- 
milab experiment E772, where X 2 is the momentum frac¬ 
tion of the parton inside the nucleus and the invariant 
mass of the produced muon pair, M, covers the range 
~ (4.5,13) GeV (excluding the charmonium and botto- 
nium resonances). These data have been taken for large 
Feynman xp ^ xi— X 2 corresponding to smallish X 2 val¬ 
ues. 

Similarly, in Fig.[^(b), we present a comparison of our 
predictions with large xp data from the E866 experiment 
for the ratio R = {da^/dxidM)/{da^/dxidM). The 
data are arranged in four bins of the invariant mass (M = 
{4.5, 5.5,6.5, 7.5} GeV) and are presented as a function 
of the proton momentum fraction xi. 

As can be seen, the theory predictions describe the 
data quite well, except for some isolated points (generally 
those with large error bars). 


3. Pion production data sets 


The newest addition to the current analysis as com¬ 
pared to Ref. [13] are the ratios of double differential 
cross-sections for single inclusive pion data from the 
STAR and PHENIX experiments at RHIC. Specifically, 
we fit the ratio 


^d^aj^VdpTdy 

dAu d'^a^^/dpTdy 


(4.1) 
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X X 


Figure 10: Comparison of the nCTEQlB NLO theory predictions for R = F^{x,Q^)/F^{x,Q'^) as a function of x 
with nuclear target data. The theory predictions have been calculated at the values of the corresponding data 
points. The bands show the uncertainty from the nuclear PDFs. 



Figure 11: Same as in Fig. 


10 


for R 


Fi{x,Q^)/Fi'{x,Q^). 
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Figure 12: Comparison of the nCTEQlB NLO theory predictions for R = jF^ as a function of with nuclear 
target data from the NMC collaboration. The bands show the uncertainty from the nuclear PDFs. 




X X 

Figure 13: Ratio of the F 2 structure functions for iron and deuteron calculated with the nCTEQlB fit at 
(a) = 5 GeV^ and (b) = 20 GeV^. This is compared with the fitted data from SLAC-E049 [57] 

SLAC-E139 ISH SLAC-E140 [59] BCDMS-85 [56] BCDMS-87 [60] experiments and results from EPS09 and HKN07. 

(The data points shown are within 50% of the nominal value.) 


and we include only the data measured at central rapid- and 0.962 for PHENIX and STAR, respectively. These 
ity to exclude potential final-state effects (this criterion 
excludes any data from BRAHMS). Additionally, we fit 
the normalizations of the RHIC data and obtain 1.031 
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Figure 14: Comparison of the nCTEQlS NLO theory predictions for R = ctqy/'^dy with data for several nuclear 
targets from the Fermilab experiments E772 (left) and E866 (right). The error bands show the uncertainty from the 

nuclear PDFs. 


values are within the experimental uncert aintyj^ Fitting 
the single inclusive pion production has the added com¬ 
plication that it depends on the fragmentation functions 
(FFs). As mentioned in Sec. [ITj pre-computed grids of 
convolutions with the free deuterium PDFs and a set of 
FFs are used to speed up the NLO calculation. 

In Fig. |15a[ PHENIX and STAR data are com¬ 
pared with predictions from the nCTEQlS fit using the 
Binnewies-Kniehl-Kramer (BKK) fragmentation func¬ 
tions [84] . As the PHENIX data are more precise than 
the STAR data, the former will have a correspondingly 
larger impact on the resulting fit. 

The EPS09 analysis [12] also used this data and we 
compare with their result in Fig. |15b| Our central pre¬ 
diction for RJyu differs from EPS09 but lies within their 
uncertainty band; however, our estimate of the PDF un¬ 
certainties differs substantially from EPS09j^ The main 
reason for this difference is the fact that EPS09 chooses to 
include the single inclusive pion data with a large weight 
(x20) to enhance its importance, and this choice leads 
to the suppression of the corresponding uncertainties. 


Another source of difference can arise from the choice 
of the fragmentation functions. The EPS09 analysis uses 
the Kniehl-Kramer-Potter (KKP) fragmentation func¬ 
tions [5S] whereas the nCTEQlS fit is based on the BKK 
FFs. To investigate the effect of different fragmentation 
functions, we have calculated R^Au using the KKP FFs 
but still using the nCTEQlB nPDFs obtained employing 
the BKK FFs (see Fig. I6a). As can be seen, the choice 
of different fragmentation functions yields only minor dif¬ 
ferences. 

In a second step, we have also performed a complete 
reanalysis of the nuclear PDFs using the KKP fragmen¬ 
tation functions in both the fit and also for the calcula¬ 
tion of RJau ^ud this is shown in Fig. I6b The use of 
the KKP FFs does not change the central prediction for 
-^dAu dut slightly changes the nPDF uncertainties in the 
high-pT region. 

In summary the use of two different sets of fragmen¬ 
tation functions, BKK and KKP, has only a minor effect 
on the resulting nPDFs. This does not exclude a possi¬ 
bility that a larger effect on nPDFs is possible if other 
fragmentation functions are used |86j . 


We note that the EPS09 analysis obtained similar normaliza¬ 
tions. 

The EPS09 analysis uses a different asymmetric definition of un¬ 
certainties given by 

(AX+)2 = ^ [max {x(S+) - X(Sl),X{S-) - XfSg), o}] " , 

k 

[XX-f = ^ [max {x(S°) - X(S+), X(5“) - X(5,-), o}] " . 

k 

To make this comparison consistent, we adopt the same definition 
when comparing with the EPS09 prediction. 


C. Fit without inclusive pion data (nCTEQ15-np) 

To further analyze the impact of the newly added in¬ 
clusive pion data and because the pion data introduce 
an unwanted dependence on fragmentation functions, we 
performed an alternative analysis which does not include 
the RHIC inclusive pion data (nCTEQ15-np). 

In Fig. [T^ we compare the results of the nCTEQlS fit 
with the ones of the alternative analysis nCTEQ15-np. 
When examining the nuclear correction factors (left pan- 
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Pt [GeV] 



Pt [GeV] 


(a) Comparison of the nCTEQlS fit with the data. The error 
bands are computed by adding the uncertainties in quadrature. 


(b) Comparison of the nCTEQlS and EPS09 fits with the data. 
The nCTEQlS error bands are computed using asymmetric 
uncertainties (MAX) to match EPS09. 


Figure 15: We display the comparison of the nCTEQlS and EPS09 fits with the PHENIX and STAR data for 
the ratio RJau- The plotted PHENIX and STAR data are shifted by our fitted normalization. 
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(a) Comparison of the nCTEQlS fit using the default BKK (blue) (b) Same as previous figure, but with a full re-analysis using the 
and the KKP fragmentation (violet) functions for the calculation BKK (blue) and the KKP fragmentation (violet) functions 

of throughout the fitting procedure. 


Eigure 16: We compare the impact of different fragmentation functions on the observable RJau- The nCTEQlB error 
bands are computed using asymmetric uncertainties to match EPS09. 


els) we see the pion data have an impact on the gluon 
PDE and to a lesser extent on the valence and sea quark 
distributions. For the central prediction, the inclusion of 
the pion data decreases the lead gluon PDF at large x 
and increases it for smaller x] the two gluon distributions 
cross each other at a; ~ 0.08. Throughout most of the 
x-range the error bands are reduced with the exception 
of X ~ 0.1 (and very small x values) where they stay 
more or less unchanged. This is precisely the range that 
is sensitive to the DIS Sn/C (and DY) data. For most of 


the other PDF flavors, the change in the central value is 
minimal (except for a few cases at high-x where the mag¬ 
nitude of the PDFs are small). For these other PDFs, the 
inclusion of the pion data generally decreases the size of 
the error band. 

In Fig. the predictions of the nCTEQlB and 
nCTEQlB-np fits are compared to the RHIC pion produc¬ 
tion data. The effect of the pion data is to increase 
for small pt and decrease it at larger px by up to 5%. 
The two central predictions cross each other at pr 4 
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Figure 17: Comparison of the nCTEQlS fit (blue) with the nCTEQ15-np fit without pion data (gray). On the left we 
show nuclear modification factors defined as ratios of proton PDFs bound in lead to the corresponding free proton 
PDFs, and on the right we show the actual bound proton PDFs for lead. In both cases scale is equal to Q = 1.3 GeV. 


GeV. This can be connected to the crossing of the gluon 
distributions in Fig. 17 (at x ~ 0.08) which is in line with 
the kinematic mapping in Fig. 


D. Constraining the PDF flavors with data 

Global analyses of PDFs necessarily include data from 
a wide variety of experiments which are differently sen¬ 
sitive to various PDF flavors. Examining the leading or¬ 
der expressions for DIS, DY, or 7r-production provides a 
simple estimate of which observable can constrain which 
PDF flavor combination. Additionally, we have to take 
into account the number of data points and their statis¬ 
tical and systematic uncertainties. All of these factors 
contribute to the x^-function; hence, we start with this 
measure to evaluate the impact of different experiments 
upon the PDF flavors. 


1. vs. the gluon parameters 


In Fig. 

and the contributions from individual experiments to this 
change as a function of the shift of selected gluon param¬ 
eters {cf i, CQ 2 } from the respective best fit val¬ 
ues. Recall that the parameters {cf c\.^} control the 
shape of the gluon PDF whereas {Cg 2 } controls the A- 
dependence of the normalization. The remaining gluon 
parameters behave in a similar manner as cf ^ and c® ^ • 
One feature that is immediately apparent is that the 
nCTEQlB minimum is not necessarily a minimum for all 
the experiments individually. For example, we see that 
the PHENIX experiment would prefer to shift cf ^ to 
larger values (^0.002) while some of the DIS experiments 
(e.g., ID=5116, NMG-96 Pb/G) prefer a lower value for 
cf I (~-0.002). Therefore, the obtained fit is a compro¬ 
mise that depends on the relative weight of the vari¬ 
ous data sets. This observation is part of the reason 


19 we compare the change of the global x^ 
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figures highlight some of the tensions between the indi¬ 
vidual data sets that the global fit must accommodate. 

On top of that, Fig. shows which experiments are 
most sensitive to the change of the underlying gluon pa¬ 
rameters. In turn the same experiments are the ones 
which have the largest impact when constraining the 
gluon PDF. Perhaps in contrast with expectations, the 
parameters analyzed in Fig.j^are mostly constrained by 
the NMC Sn/C data and data from several other DIS ex¬ 
periments. We also see that the inclusive pion production 
from PHENIX is sensitive to the gluon shape parameters 
(cf i, c® ]^) but not to its normalization (cg 2 ). 


2. Correlations between data sets and PDFs 


Figure 18: Comparison of the predictions of the 
nCTEQlB (solid blue) and nCTEQ15-np (dashed gray) fits 
to inclusive pion production data from PHENIX and 
STAR demonstrating the effect of including these data 
sets. Note that, the dark blue area is the overlap 
between the blue and gray bands. 


we consider a = 1 tolerance crit erion impractical 
and choose Ay^ = 35 (see Appendix Al). Moreover, 
for some experiments there may not even be a local min¬ 
imum in the vicinity of the nCTEQlB solution. Thus, these 


Looking at the dependence of the y^-function on only 
three gluon PDF parameters cannot give a complete pic¬ 
ture and neither would inspecting the behavior for all 
gluon parameters because the momentum sum rule con¬ 
nect in fact all PDF flavors together. Therefore in the 
following we use different methods to study the impact 
of individual experiments on different PDF flavors. 

We introduce two quantities which will help us analyze 
the impact individual experiments have on constraining 
given PDF flavors. The first quantity is the cosine of 
the correlation angle between two observables X and Y 
which was used in 01117] and can be defined as 


cos 4’[X,Y\ 



(4.2) 


where the indices ipdf run over the 16 eigenvector directions. 

In the following we will use the cosine of the correlation angle to investigate the correlations between the y^ 
functions of the individual experiments and a single PDF. For example, in the case of the gluon PDF the cosine of 
the correlation angle has the form cos(j)[g{x,Q),x^{jexp)]- This correlation cosine depends on x and Q through the 
gluon PDF, g{x,Q), and on the particular experiment through X^Uexp)- 

Even though the cosine of the correlation angle is a useful quantity, it doesn’t highlight the experiments with more 
data or smaller errors. It turns out that the normalization factors in Eq. (4.21 strongly reduce any sensitivity to the 
number of data points or to the size of the errors of an experimental data set. Therefore we introduce an alternate 
measure, the effective y^ for an experiment jexp, defined as 


^Xesidexp, X) — 


^ 2 

ipdf 


X: 


2 (+) 
ipdf 


Uexp) Xipj^f (jexp) 


,r Uexp) Xi 


^‘^pdf 


2 ( 0 ) 

ipdf 




(+1 — x ^~'^ 

ipdf ipdf 



As before, the index ipdf runs over the 16 Zi eigenvector 
directions. 

Ay^fj is positive definite and comparing the definitions 


(4.2) and (4.3) it is missing the normalization factor for 
the y^ function which allows it to be more sensitive to 
experiments with more data or smaller errors, i.e., exper- 










































































23 



total y 

» * DY; W/Be (FNAL-E886-99) 

■ ff" : DAu/pp (STAR) 

•—* DIS; Pb/C (NMC-96) 

* DIS; Sn/C (NMC-96) 

Other Experiments 

■ / : DAu/pp (PHENIX) 




cl^ ( 0 . 05488 ) 


total y 

• « DY: W/Be (FNAL-E886-99) 

■—■ tt” : DAu/pp (STAR) 

« DIS: Pb/C (NMC-96) 

* DIS; Sn/C (NMC-96) 

Other Experiments 

: DAu/pp (PHENIX) 



100 



-0.015 -0.010 -0.005 

C® 

0, 

0.000 0.005 0.010 0.015 

(-0.0373) 

total y^ 

« « DY: W/Be (FNAL-E886-99) 

• DIS: Ca/D (NMC-95,re) 

□ □ DY; Fc/H2 (FNAL-E772-90) 

i--* DIS; Sn/C (NMC-96) 

DY; W/H2 (FNAL-E772-90) 

» DIS; C/D (NMC-95) 

Other Experiment.? 


Figure 19: Contribution of different experiments to the total ~ Xo function (solid black line) for a 

selection of gluon parameters (a) cf (b) cf (c) Cg 2 - On the cc-axis we show the shift from the best fit value 

(indicated in the parenthesis, cf. Table [v|). 


iments which have a larger impact in constraining single 
PDF flavors. 

we display both the and cor¬ 


In Figs. 20 and 
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relation cosine as a function of x. These plots do not 
exhibit a strong Q dependence, so we only display them 
for one value of Q = 10 GeV. 

We now examine the ^XeS results for the gluon PDF 
in lead (A=207) as show in Fig. 20a[ For readability, 
we primarily show the data sets which have the largest 
impact on Ay^g; these are generally the data sets which 
involve the heaviest targets. The strong influence of the 
DIS Sn/C set reflects a combination of the large cov¬ 
erage of the data and the small errors. The DIS Pb/C 
data, and to a lesser extent the DIS Sn/D data, also pro¬ 
vide constraints for the gluon PDF in lead. The PHENIX 
pion production data contributes strongly in the central x 
region; conversely, the effect of the STAR data is negligi¬ 
ble due to the larger uncertainties. Additionally, the DY 
data on heavy targets (W tungsten with Be and D) also 
play a role in determining the gluon lead PDF; this is due 
to the fact that the DY data cover a range ^ (20,170) 
GeV^ in the invariant mass of the muon pair, which cre¬ 
ates some sensitivity to the gluon PDF via scale evolu¬ 
tion. 

In Fig. |20b| we show the correlation cosine for the gluon 
PDF in lead. The DIS Sn/G and DY W/Be data sets 
have positive correlations at large and small x, and a neg¬ 
ative dip in the middle. Gontrary to this, the DIS Pb/G, 
Sn/D and DY W/D data sets have the opposite behavior. 
Hence, these data sets are anti-correlated which indicates 
that they pull against each other in the fit. This is pre¬ 


cisely what we have observed in Fig. 19 for the gluon 
parameters. Also, the PHENIX data have a separate 
cc-dependence (arising from a separate production mech¬ 
anism), and this will further help us separate the PDF 
flavor components. 

Finally, there are two data sets (STAR and DIS Xe/D) 
that have relatively large correlation cosines, but do not 


have a large influence on the Ay^g of Fig. 20a thus, we 
need to take care when interpreting the results of the 
correlation cosine plots and use this in combination with 

We now consider the gluon PDF in carbon (A=12) to 
see if the general observations above apply in the case of 
a lighter nuclei. In Fig. |20c| we see the primary data sets 
constraining Ay^g are the DIS sets involving ratios of 
carbon (Sn/G, G/D, Pb/G) or other comparable nuclei 
(Ga/D). Note the DY data on heavy tungsten (W) and 
the pion production data on gold (Au) are not shown as 
they do not contribute significantly to Ay^g for carbon. 

The correlation cosines for the gluon PDF in carbon 
are shown in Figs. 20d We see the DIS Pb/G data have a 
positive correlation cosine at small x and a negative cor¬ 
relation cosine at large x. The DIS Sn/G data shows the 
opposite behavior; hence, these data sets will pull against 
each other in the fit. The DIS G/D and Ga/D data gen¬ 
erally have a small correlation cosine throughout the x 
range. As in the case of the gluon in lead, we see there 
are a number of data sets (such as DIS Fe/D) that have a 
large correlation cosine but yield a small contribution to 
the Ay^g; thus, we need to use both the Ay^g and cos cj) 
information together when drawing our conclusions. 

We now turn our attention to and distri¬ 

butions for lead at Q = 10 GeV as shown in Fig. [21] 
For these PDFs, not only is the Q-dependence rather 
mild, but the differences between heavy and light nuclei 
are also not as pronounced as in the gluon case. The 
Ay^g for the u and d PDFs depends almost exclusively 
on the DIS data from heavy targets (Sn/G, Pb/G, Fe/D), 
with some contributions from PHENIX pion production 
data at small-cc, and a minimal contribution from the DY 
W/D data. 

Turning to the coscj) plots, we see the DIS Sn/G and 
the Fe/D data start with a positive correlation cosine 
at small x and moves negative for increasing cc, while the 
DIS Pb/G and the DY W/D data do the opposite; hence. 
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Other Experiments 


(a) gluon: Ax^f^{g) at Q = 10 GeV for lead 




(b) gluon: cos(l>{g,x^) at Q = 10 GeV for lead 





(d) gluon: cos0(p,x^) at Q = 10 GeV for carbon 


Figure 20: Correlation measures for lead and carbon at Q = 10 GeV for the gluon of the nCTEQlb fit. The left 
panels display the effective and the right panels display the correlation cosine as a function of x. 


these sets are anti-correlated in this region. For small to 
medium x values, the general pattern is similar between 
the u and d correlation plots, but they differ some at 
large x where we see, for example, the DIS Fe/D data 
has a positive correlation cosine for u but a negative one 
for d; this will be useful in differentiating u and d PDFs 
at large x. As with the gluon correlation plots, there are 
a number of data sets (such as the DIS Ag/D) which have 
large correlation cosines but small contributions to 
thus, they have minimal effect constraining the PDFs. 


E. Comparison with different global analyses 

We now compare our nCTEQlS PDFs with other recent 
nuclear parton distributions in the literature. Specif¬ 


ically, we will consider DSSZ [TT], EPS09 [T^, and 
HKN07 Our data set selection and technical as¬ 

pects of our analysis are closest to that of EPS09. In 
Eigs. and we plot nuclear modifications for the 
PDEs of a proton bound in lead, pIP^jfv (left), as well 
as the bound proton PDFs themselves, (right), for 

different flavors for a selection of Q scales. 

For the u and d PDFs at Q = 2 GeV, nCTEQlS has 
significant overlap with the other sets through much of 
the X range with a stronger shadowing at small x. Our 
results at X < 10“^ are extrapolated since they are not 
constrained by data due to the cut Q > 2 GeV which was 


Note that there is also a very recent global nPDF analysis per¬ 
formed at NNLO level 1881 . 
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(a) ti-quark: Ax^^{u) at Q = 10 GeV for lead 
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(b) li-quark: cos4>{u,x^) at Q = 10 GeV for lead 
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(c) d-quark; Ax^^{d) at Q = 10 GeV for lead 
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(d) d-quark: cos0(d,x^) at Q = 10 GeV for lead 


Figure 21: Correlation measures for lead at Q = 10 GeV for the ii-quark and d-quark distributions of the nCTEQlS 
fit. The left panels display the effective and the right panels display the correlation cosine as a function of x. 


imposed in order to reduce higher twist contributions. 
Therefore, it is likely that the uncertainty band at a: < 
10“^ underestimates the true PDF uncertainties. While 
this trend repeats itself for the strange quark PDF, the 
spread at small x is slightly increased!^ In fact, at Q = 
2 GeV the small x behavior of the strange PDF of all 
four fits is quite distinct with little overlap between the 


In this analysis the s-quark nuclear effects are completely deter¬ 
mined by the u and d nuclear PDFs and by the gluon nuclear 
PDF through evolution. Due to these constraints the error of 
the s-quark nuclear PDF is underestimated. A comprehensive 
analysis would require including the charged-current i^-DIS data 
as in m along with using a proton PDF baseline where the 
strange distribution was determined from different data such as 
the W c production at the LHC. 


uncertainty bands (see Fig. [2^). As we move to higher Q 
values, the DGLAP evolution tends to bring the various 
PDF sets into closer agreement, particularly at small x 
values. For example, already at Q = 10 GeV the nCTEQlS 
bands overlap the other PDFs across a much broader x 
range than at low Q values. 

In the case of the gluon, there is considerable variation 
among the different PDF sets at Q close to the initial 
scale. Again, the nCTEQlB exhibits a stronger shadowing 
suppression along with a larger enhancement in the anti¬ 
shadowing region (x ~ 0.1). In addition, the uncertainty 
band for x > 0.02 is considerably larger than the uncer¬ 
tainty bands of the other groups. The nCTEQlB result is 
largely compatible with the result of EPS09 even though 
the shape of the central prediction is more suppressed in 
the shadowing region and enhanced in the anti-shadowing 
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Q=2 GeV 


Figure 22: Comparison of the nCTEQlS fit (blue) with results from other groups: EPS09 [12] (green), DSSZ [TT] 
(orange), HKN07 [14] (red). The left panel shows nuclear modification factors for lead, and the right panel the 
actual PDFs of a proton bound in lead . The scale is Q = 2 GeV. The wide spread of the ratios at large x are an 
unphysical artifact due to the vanishing of the PDFs in this region. 


region. We have less overlap with the HKN07 and DSSZ 
bands, in part, due to their smaller uncertainty bands. 
Moving to larger Q values, the DGLAP evolution again 
causes the different PDFs to converge,. 

Note that the ratio plots of Figs. [2^ and [^ have quite 
a wide spread at large x values. This unphysical behavior 
is an artifact due to the vanishing of the PDFs in this re¬ 
gion. The spread is largest for those PDFs with minimal 
support at large a;-specifically g, s, u, d. Also, these ef¬ 
fects are reduced when we construct the full nuclear lead 
distribution as shown in Fig. [24] 

Examining the u- and d-valence distributions, one can 
see that PDE sets {HKN07, EPS09, DSSZ} agree quite 
closely with each other throughout the x range. While 
the nCTEQlB fit uncertainty bands generally overlap the 
other sets, we see on average the Uy distribution is softer 
and the dy distribution is harder. These differences re¬ 
flect the fact that the HKN07, EPS09, and DSSZ fits 


assume that the nuclear corrections Ru^ and Rd^ are the 
same, while the nCTEQlB fit allows them to vary inde¬ 
pendently. Clearly, there is no physical reason to assume 
that Uy and dy must have a universal nuclear correction 
factor, and there exist models in the literature [zniiHniiHi] 
which indeed predict non-universal modifications. 

The obvious question is whether the additional free¬ 
dom to decouple the and Rd^ nuclear corrections 
yields a substantial improvement in the fit. To shed more 
light on this issue, we have generated a modified fit where 
we have forced the Uy and dy nuclear corrections to be 
similar to the EPS09 PDE setj^ We find that the y^/dof 
for this modified fit is comparable (Ay^ < 5) to our orig- 


As we are fitting directly the nuclear PDFs and 

not the ratios f^^^{x,Q)/f^{x,Q), it is non-trivial to force 
the nuclear corrections to be exactly the same if the under- 
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Figure 23: Same as Fig. 


22 


with Q = 10 GeV. 


inal nCTEQlB at a level well below our tolerance criteria 
of = 35. Therefore, we conclude that the current 
data sets are not sufficiently sensitive to distinguish the 
Ui, and dy nuclear corrections to a good degree. Hence, 
the advantage of independent and Rd^ correction 
factors is currently limited, which however will change 
with more data (e.g. from the LHC)(f] 

To better understand this result, we observe in Figs.[22| 
and that the and dy ratios exhibit opposite 


lying proton PDFs differ. We are able to find an approxi¬ 
mate solution by equating the and dv coefficients Cij for 
{ij} = {11,12, 21, 22, 31, 32, 51, 52} and refitting the PDFs’. 

In an earlier study we did find an apparent difference due to 
independent Ru^ and Rd^ nuclear corrections. The present up¬ 
dated analysis additionally includes; i) an improved treatment 
of the {A, Z} isoscalar corrections, ii) QED radiative corrections 
for DIS data sets, iii) use of full theory (instead of K-factors) to 
obtain the final minimum, and iv) improved numerical precision 
for the DY process. With these improvements, the 
modified fit is now comparable to nCTEQlS. 


x-dependence as compared with the {HKN07, EPS09, 
DSSZ} sets. That is the Uy ratio is below the other sets 
at large x and above at small x\ the dy ratio does the 
opposite. As the nuclear data sets probe a linear combi¬ 
nation of Uy and dy , this raises the question as to whether 
the above differences might cancel when combined. 

In Fig. [^we now compare the full nuclear lead PDFs 
from the different groups. The upper panel shows the 
PDFs themselves, and the lower one shows their ratio 
compared to the nuclear combination constructed out of 
the free proton - the full nuclear correction. From this 
comparison we can clearly see that the large differences in 
the effective bound proton distributions of valence quarks 
(Figs. 22 231 translate into much smaller differences in 
the full nuclear PDFs that actually enter the calculation 
of observablesIn particular, we see that Uy and dy 


Note that only up and down distributions differ between the full 
nuclear PDFs and the PDFs of the bound proton; the gluon and 
strange distributions are the same. 
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Figure 24: (upper panel) Comparison of the full nuclear lead distributions, + ^° 207 ^^ jn/Pb^ 

nCTEQlB (blue), EPS09 (green) and HKN07 (red) at Q = 10 GeV. Lower panel shows the same distributions 
compared to the lead PDF, constructed of free proton distributions. The wide spread of the ratios at large x 
are an unphysical artifact due to the vanishing of the PDFs in this region. 


distributions of the nCTEQlB fit are in very good agree¬ 
ment with the EPS09 results, and have substantial (but 
not complete) overlap with HKN07p^ 

Of course, as the data can only constrain the full nu¬ 
clear PDF in the combination = \pl^ + P^^i 

we conclude that better separation of Uy and dy distri¬ 
butions require more data on non-isoscalar targets. We 
also note that the currently available DIS data use a num¬ 
ber of non-isoscalar targets and would have the potential 
to partially distinguish Uy and dy distributions; unfortu¬ 
nately many of these data sets have been corrected for 
the neutron excess and in turn lost this ability. 


The DSSZ set (not show) is similar to HKN07 in that it has 
substantial (but not complete) overlap. 


V. SUMMARY AND CONCLUSIONS 

In this paper we have presented the first complete anal¬ 
ysis of nuclear PDFs with errors in the CTEQ framework. 
The resulting fit, nCTEQlB, uses the available charged lep¬ 
ton DIS, DY and inclusive pion data taken on a variety of 
nuclear targets. The uncertainty of this analysis is pre¬ 
sented in the form of error PDFs which are constructed 
using an adapted Hessian method. 

Within our framework we are able to obtain a good 
fit to all data. The output of the nCTEQlB analysis is a 
complete set of nuclear PDFs with uncertainties for any 
A = {I,..., 208}. A selection of nuclear PDFs for the 
most common nuclei are made publicly available but 


The nPDF sets for the current nCTEQlS analysis as well as for 
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custom nPDFs can be generated for any {A, Z} combi¬ 
nation. 

In comparison to our previous analysis m, we have 
included the data from the inclusive pion production 
at RHIC. The new data provide additional constraints 
mostly for the nuclear gluon PDF but the description of 
the data relies on the fragmentation functions. There¬ 
fore we also provide an alternative conservative result 
nCTEQ15-np which does not include the inclusive pion 
data and is hence fragmentation function independent. 

Compared to other global analyses (HKN07, EPS09, 
and DSSZ) there are a number of important differences: 

• In contrast to the other analyses, we parameterize 
the nuclear PDFs directly instead of the nuclear 
corrections factors. 

• In addition, our u- and d-valence distributions are 
parametrized independently. 

• Other differences arise from the selection of data 
points used in the fit. In particular we impose more 
conservative kinematic cuts in order to minimize 
effects from higher twists and target mass correc¬ 
tions. 


Overall our results are compatible with the other 
nPDFs but after a detailed look we see distinct differ¬ 
ences (see Fig. 24). 


(i) The nCTEQlB nuclear gluon PDF has a larger shad¬ 
owing at small-x than the other global analyses. 
Our result is compatible with the result of EPS09 
as the error bands are overlapping throughout the 
entire x range. The overlap in case of HKN07 and 
(especially) DSSZ is limited especially in the small¬ 
er region where no data constraints are present (and 
uncertainties of HKN07 and DSSZ are very small). 
This highlights the fact that nPDF uncertainties, 
in particular for gluon, are underestimated and dif¬ 
ferent gluon solutions are possible [33]. 


(ii) Our valence distributions for a bound proton in 
lead differ as we allow separate nuclear corrections 
for Uy and dy. Compared to the other groups; 
our d-valence PDE is harder and u-valence PDF 
is softer. However, when the full lead nucleus is 
constructed, these differences are substantially re¬ 
duced and we observe a good agreement between 
all groups. 


(hi) The nCTEQlS light sea quark distributions are in 
very good agreement with the ones from the other 
groups for x > 10“^. At smaller x where there 
are no data constraints the individual error bands 
clearly underestimates the uncertainty. 


the alternative nCTEQlS-np analysis are available for download at 
http://ncteq.hepforge.org as well as on the LHAPDF website. 


(iv) It should be also mentioned that strange distribu¬ 
tions are currently not fitted in any of the nPDF 
analyses and are fixed by imposing additional as¬ 
sumptions; this leads to quite significant differences 
between different groups. 

All in all we find relatively good agreement between 
different nPDFs. Most of the noticeable differences be¬ 
tween them occur in regions without any constraints from 
data and so they can be attributed to different assump¬ 
tions such as parameterization of the nuclear effects. 

In view of the differences, the true nPDF uncertainties 
should be obtained by combining the results of all analy¬ 
ses and their uncertainties. In particular, this is true for 
the gluon distribution where the small x behavior is basi¬ 
cally unconstrained and every single nPDF analysis sub¬ 
stantially underestimates it (see our earlier study [33]). 

The nCTEQ framework used for the nCTEQlS fit can 
combine data from both proton and nuclear targets into 
a single coherent analysis. Using nCTEQlB fit as a refer¬ 
ence, it will be interesting to include the upcoming LHC 
data as we continue to investigate the relations between 
the proton and the nuclear PDFs. 
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Appendix A: Determination of Ay^ and Hessian 
rescaling 

1. Determination of Ay^ 

In this appendix we discuss the details of the determi¬ 
nation of Ay^ which is motivated by the treatment pre¬ 
sented in Refs. [El SSI ST). We investigate how the global 
fit describes each experiment by examining y^ which is 
the individual y^-contribution of experiment k with Nk 
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Figure 25: The 90% confidence level limits from different data sets in the eigenvector direction zi. The x^-minimum 
for each experiment is denoted by a black square, and the green band demonstrates the interval of the eigenvector 

parameter corresponding to the final 


data points. We can then see how xl changes when vary¬ 
ing PDF parameters along each eigenvector direction Zi 
of Eq. (2.25). 


The probability distribution for the Xk given that the 
fit has degrees of freedom is: 


PixlNk) = 


(^2)W./2-lg-x=,/2 

2Nk/2Y(Nk/2) 


This allows us to define the percentiles via 


(Al) 


/ Pix'^j^)dx^ = P% where p = {50, 90, 99} . 

Jo 

(A2) 

Here, ^50 serves as an estimate of the mean of the dis¬ 
tribution and ^ 90 , for example, gives us the value where 
there is only a 10 % probability that a fit with x^ > Cqo 
genuinely describes the given set of data. 

Due to fluctuations in the data and possible incom¬ 
patibilities between experiments, the global x^ minimum 
does not necessarily coincide with y^-minima of individ¬ 
ual experiments. Moreover, for the same reason, the min¬ 
imum x^ for each experiment, Xk qj f®"'^ away from 

the expected minimum given by ^ 59 . In order to use the 
percentiles defined in Eq. (A2) to dehne the 90% confi¬ 


dence level (C.L.), we rescale the ^go percentile to take 
into account the position of the minimum as 


where the Xk stays within the 90% C.L. limit (i.e. Xk < 
^go). Eor each eigenvector direction we then construct an 
interval {z~,zJ) where all experiments stay within the 
90% C.L. limit as 

(A5) 

k 


These intervals can obviously be different for each eigen¬ 
vector, depending on the fact how well the experi¬ 
ments constrain the variations in this eigenvector direc¬ 
tion. Eor n free parameters we obtain 2n parameters 
{z^, z^, Z 2 , Z 2 , ■ ■ ■, z~, z^} which we can use to define 
the global tolerance as 


i 


2n 


(A 6 ) 


Having performed the procedure described in this Sec¬ 
tion, we have arrived at Ay^ = 35. One can compare 
how this choice of global tolerance (the same for every 
eigenvector direction) agrees with the rescaled 90% confi¬ 
dence level (C.L.) for each experiment in every direction. 
In Fig. we show this comparison for only one single 
eigenvector direction as all the others are rather similar. 


2. Hessian rescaling 


^90 —^ ^90 

For each eigenvector direction given by a variation of 
the parameter Zi and every experiment, we define an in¬ 
terval 


(fc)— ^ ~ ^ {k)+ 

zi <zi<zy , 


(A4) 


Choosing a larger tolerance Ay^ = 35 as argued above 
might pose a problem for the Hessian approach as it 
requires using information from a larger neighborhood 
of the global minimum which is not necessarily well de¬ 
scribed in the quadratic approximation. Fig. |26| confirms 
that this is the case for the nCTEQlB fit. Both in the 
original pa rame ter space, Fig. |26a| and in the eigenvector 
basis. Fig. 26b we can see directions where y^-function 
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(a) function in parameter space. 


(b) x^ function in eigenvector space. 


Figure 26: These plots display the Hessian before the “rescaling” procedure, 
function relative to its value at the minimum, Ax^ = ~ Xo^ plotted along the 16 fitting parameters of the 

original space (left) and along the zt directions in the eigenvector space (right). The actual function is plotted 
with solid lines, and the Hessian approximation Ax^ = zf is shown with dashed lines. 


deviates substantially from the quadratic approximation 
when Ax^ ^ 35. This is a problem because in the Hes¬ 
sian approach we use the eigenvector basis to determine 
the ranges of the normalized parameters Zi Fig. |26| shows 
that if we take Azi = Ax^ « then depending 

on the specific eigen-direction we would largely overesti¬ 
mate or underestimate the error on our parameters (see 
e.g. plots #1, #2 and #14 in Fig. 26b). 

To improve the constraints provided by the x^- 
function, we redefine the Hessian which we use to de¬ 
termine the error PDFs using the formalism described in 
Sec. |HC| We keep the eigenvector information intact, but 
rescale the eigenvalues of the original Hessian (which cor¬ 
responds to rescaling the parameters Zi) so that the mod¬ 
ified Hessian better describes the y^-function not only in 
the minimum (Ay^ = 0) but also at = 35. For each 
eigenvector direction, we identify the parameter values 
zf where 


^X^{zf) = X^{zf)-xl = i^, (A7) 

where Xo is the minimum of the x^- Using the z^, we 
rescale the corresponding eigenvalue as 


Ai I—> A( — 


15+|U 


?-|2 


2yAp 


A, 


(AS) 


The impact of the rescaling of the Hessian can be seen 


on Fig. 27b where one notices that the description of the 


X^-function in the eigenvector basis is highly improved, 
especially in region where Ay^ = 35. The descripti on of 
the x^-function in the original parameter space (Fig. 27a) 
is also improved but to a lesser extent. However, this is 
a secondary feature as we are working in the eigenvector 
space when defining the error PDFs. 


Appendix B: Usage of nCTEQ PDFs 

We provide a set of PDF tables for the nCTEQlS and 
nCTEQ15-np hts at the nCTEQ Hepforge website [HS]. We 
provide the tables in the older CTEQ PDS format to¬ 
gether with a dedicated interface as well as in the new 
LHAPDF6 format [20] ■ In the future the LHAPDF6 
grids will be also available at the LHAPDF website m- 

We provide tables for both bound proton PDFs 
fp/^[x,Q) as well as grids for the resulting full nuclear 
PDFs = Z/A + {A - Z)IA /"/^. The bound 
proton PDFs allow a direct comparison of the nPDFs for 
different A values as displayed in Fig. On the other 
hand the full nuclear PDFs can be used directly to cal¬ 
culate cross-sections for the nuclear collisions. 

At the moment we distribute grids for a selection of nu¬ 
clei that are commonly used in the high energy/nuclear 
experiments. In particular, we provide girds for: {He, 
Li, Be, C, N, Al, Ca, Fe, Cu, Kr, Ag, Sn, Xe, W, Au, 
Pb}. Since our parametrization is continuous in A and Z 
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(b) function in eigenvector space. 


Figure 27: These plots display the Hessian after the “rescaling” procedure, 
function relative to its value at the minimum, ~ Xo^ plotted along the 16 fitting parameters of the 

original space (left) and along the Zi directions in the eigenvector space (right). The actual function is plotted 
with solid lines, and the Hessian approximation = zf is shown with dashed lines. 


it allows us to generate PDFs for any nuclei or isotopes. 
In case users are interested in having the nCTEQlS distri¬ 


butions for a nucleus that is not included in our standard 
selection, we can generate the PDFs upon request. 
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